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Abstract 

A new and efficient algorithm is presented for the calculation of the 
partition function in the 5 = ±1 Ising model. As an example, we use 
the algorithm to obtain the thermal dependence of the magnetic spin 
susceptibility of an Ising antiferromagnet for a 8 x 8 square lattice 
with open boundary conditions. The results agree qualitatively with 
the prediction of the Monte Carlo simulations and with experimental 
data and they are better than the mean field approach results. For 
the 8x8 lattice, the algorithm reduces the computation time by nine 
orders of magnitude. 

Keywords: antiferromagnets; Ising model; magnetic susceptibility; mean field 
approximation; world records. 

1 Introduction 

A calculation of the partition function exactly in the thermodynamic limit is 
equivalent to knowledge of all equilibrium properties of a given system. This 
task is usually beyond our horizon, except for some selected cases which we 
call "trivial" . The Ising model seems to be at the border of our possibilities 
since the beginning of the 20-th century. Even for a finite system of iV spins 
our computational ability is limited by the number of the system configura- 
tions, which is 2^. Here we describe a new algorithm, designed to improve 
the speed of the calculation of the partition function for a finite Ising system. 
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Table 1: The CPU time necessary for investigation of M different configu- 
rations of 8 x 8 large 2D Ising lattice on SGI 2800 machine. 

~m W W W 10 s To 5 W° 

t c pu [sec] 0.86 6.90 66.4 660 6648 65752 



In the Ising model P only two spin states are possible, say "up" and 
"down" [Si = ±1). The energy E of a given spin configuration may be 
expressed in terms of the number n of spins pointing "up" (say Si = +1) and 
the number k of anti-parallel bonds between the nearest neighbors (SiSj = 

-1): 

E(n, k) = -jJ2 S * S i ~ H J2 S * = 2J ( k ~L 2 + L )- H (L 2 ~ 2n), (1) 

<i,j> i 

where J denotes the exchange integral, H is an external magnetic field, and 
L — the linear size of the lattice. 

The partition function can be written as 

Z = ^n(n,k)-exp[-pE(n,k)], (2) 

n.k 

where Q(n, k) is the number of lattice configurations with given numbers n 
and k, 1/(3 — ksT, ks is a Boltzmann constant and T is temperature. Then, 
the average value of any quantity A may be calculated as 

(A) = Z- 1 A ( n , k ) ■ fi K k ) ■ exp[-pE(n, k)]. (3) 

n,k 

The magnetic susceptibility per spin x m ay be also expressed in the terms 
of n, k and L: 

X = /9[<5?> - (S*) 2 } = /3[((2n - L 2 ) 2 ) - (2n - L 2 ) 2 }, (4) 

where the index k enters through the averaging procedure. 

A computation of the partition function (J2J requires an evaluation of the 
histogram Q(n, k). Tab. ^shows the CPU time needed to check the numbers 
n and k of M configurations of a 8 x 8 large lattice on SGI 2800 machine. 
A rough estimation of the CPU time necessary for full investigation of all 
of 2 64 pd 10 19 possible lattice configurations gives the value larger than four 
millions years — what makes traditional/direct method practically useless. 
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Figure 1: Subsequent stages of the computation of the histogram Q$ x s(n, k) 
described in Sec. 12.11 



2 Calculations 
2.1 The algorithm 

An effective way of generation Q 8x8 (n, k) is a successive merging of smaller 
lattices and their histograms, namely f2 4x4 (n, k) and f2 8x4 (n, k). However, 
the procedure requires storing information on the Q dependence not only on 
n and k, but also on b r — the state of the r-sites-long lattice border. 

In Fig. ^a) a L 2 -long bit-string equivalent to L x L large array (see Fig. 
IHb)) is presented. For L = 4 the bit-string corresponds to an integer number 
from the interval [— 2 15 ,2 15 ). This correspondence allows to investigate all 
possible lattice configurations in a simple manner. Dark sites in Fig. E^a) 
correspond to the dark border of the lattice in Fig. Q^b), and they can be 
represented by an integer number < b 7 < 127. The first four bits of b 7 
(marked as dark sites in Fig. [He)) allow to determine the additional number 
of anti-parallel bonds, created by merging two 4x4 lattices together to get 
a 8 x 4 lattice (see Fig. ^f)). On the other hand, the last four digits of b 7 
(dark sites in Fig. HJd)) are equivalent to the part of the border < b 8 < 255 
(the dark sites in Fig. dg)) of 8 x 4 lattice. 

In the next step, two 8x4 lattices are merged to create the desired 8x8 
large lattice, as presented in Fig. [TJh). 
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Histograms (1^x4 and f2 8x8 ) for the lattices larger than 4x4 may also 
be easily computed basing on fl^x^. 

fi 8x4 (6 8 ,ni + n 2 , fa + k 2 + k') = ^2 ^4x4(&I,ni, fci) • fi 4x4 (62,n 2 , fc 2 ), 

b~2,ri2,k2 

where < k! < 4 is the additional number of anti-parallel bonds in the darker 
part of Fig. QJf) and b 8 is combined from b\ and b\. 

The histogram fls x s(n, k) may be constructed in a similar way: 

^8xs(™l + Tl 2 , fa + k 2 + k") = ^8x4(6l,7ll, fa) ■ ^8X4^2) n 2, fa), 

bf,n\,ki 
fe| ,ri2 ,A:2 

and again < k" < 8 is the number of anti-parallel bonds in the darker part 
of Fig. mh). 

This procedure on SGI 2800 machine takes only 22 hours of the machine 
time instead of a few million years in case of the usage of the traditional/ direct 
method. Successive merging may be repeated recursively to obtain the par- 
tition function for larger lattices. 

2.2 Monte Carlo simulation, mean field approach and 
experimental data 

To evaluate the results obtained for the 8x8 lattice, we show also the data 
obtained by the Monte Carlo simulations, the results of the mean field model, 
and — last but not least — the experimental data. Standard Monte Carlo 
Metropolis algorithm [2] is applied to determine the magnetic spin suscepti- 
bility of a 1000 x 1000 Ising lattice. After getting equilibrium, each point of 
the plot is obtained as the time average over a thousand of time steps. 

The mean field model is a direct generalization of the case of a ferromag- 
net. Namely, we solve numerically a set of two equations for two sublattices 
a and 7: 

m a = tanh (/5(Jm 7 + H)) 
m 7 = tanh (3(Jm a + H)) 

where J < 0. In this model, the Curie temperature Tc = —J/ks and the 
susceptibility x is found as (m a +m 7 ) / H for a small value of applied magnetic 
field, e.g. H = 0.001. 

The experimental data are collected from Ref. 0. They concern two- 
dimensional Ising antiferromagnets Rb2CoF4 and K2C0F4 where S — ±1. 
The Van Vleck susceptibility is subtracted to obtain a pure spin contribution. 
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3 Results and discussion 



In Fig. 121 the histogram Qg x g(n,k) obtained by means of the method de- 
scribed in Sec. 12. H is presented. The numbers are available at our web page 
[I]. The symbols on n-k plane (see Fig. 12(a)) indicate (n, k) pairs with non- 
zero values of f2 8x8 (n, k). In Fig. |2[b) the number of configuration f2 8x8 (n, k) 
is presented. For fixed n different symbols correspond to different k. 

In Fig. El we show the thermal dependence of the antiferromagnetic 
(J < 0) susceptibility, defined in Eq. (jlj). The Curie temperature is assumed 
to be at the inflexion point of the curve, which is at the half of the peak 
height. Below the Curie temperature, the values of x obtained with our new 
algorithm fit well the results obtained with the Monte Carlo results and the 
experimental data. Above Tc, the agreement is only quantitative. We hope 
that it can be also quantitative if periodic boundary conditions are applied. 
However, in this case the lattice border is longer, and so is the computation 
time. Still, even with the present method the qualitative accordance is better 
than the results of the mean field theory. 

It would be of interest to apply the finite size scaling to our results to 
evaluate the critical properties of the system. Then, the results could be 
compared with other methods, e.g. the finite size scaling renormalization 
group [5] or the phenomenological renormalization For this purpose, 
however, periodic boundary conditions seem to be more appropriate as a 
starting point of the calculations. 

The function Q(n, k), once known, can be easily used for the calculation of 
all equilibrium thermodynamic properties, for ferro- and antiferromagnets, 
various values of temperature and magnetic field. The summation over n 
and k is much faster, than the summation over 2 64 spin configurations. The 
results can be relevant also for other applications of the Ising model, e.g. for 
the percolation problem. As for our knowledge, the partition function has 
never been calculated exactly for the lattice as large as 8 x 8. In principle, 
the algorithm can be applied to larger lattices, with a cost of more time and 
memory. 

The computational mountain remains infinite, but its slope is a little bit 
reduced. 
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Figure 2: The histogram Q 8x8 (n, k). 
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Figure 3: The antiferromagnetic susceptibility x normalized to its maxi- 
mal value Xmax, as dependent on temperature T normalized to the Curie 
temperature T c . 
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